# Load Necessary Packages
library(glmnet)
library(foreign)
library(estimatr)
library(dplyr)
library(tidyr)
library(bcf)
library(dbarts)
library(texreg)
library(Hmisc)
library(grid)
library(gridExtra)
library(interplot)
library(margins)
library(stringr)
library(factoextra)
library(TOSTER)
library(meta)
library(dmetar)
library(interflex)

# Load Data
setwd("/users/josephphillips/Dropbox/COVID-19/Data/Canada COVID Study")
canada <- read.dta("misinfo_rd2_nosl.dta")

# Load GLMNet Function
run_LASSO <- function(data, yvar, pre){
  
  set.seed(2334)
  
  var_quo <- enquo(yvar)
  var_quo1 <- enquo(pre)
  
  ## select x variables 
  dat <- data %>%
    dplyr::select(university, agegroup, male, married, frequentchurch, region, left, right, ideology1, knowledge, nonwhite, polinterest,trust_healthgov,total_media_trust,hotspot)
  
  pre <- data %>%
    dplyr::select(!!var_quo1) %>% pull() # Save as a vector instead of a data frame ot use model matrix
  
  ## select y-variable
  out1 <- data %>% 
    dplyr::select(!!var_quo) %>% pull() # save as a vector instead of a data-frame to use model.matrix()
  
  merged1 <- cbind(out1, dat, pre) %>% na.omit() ## merge and delete NA
  x1 <- model.matrix(out1 ~ ., data = merged1) ## make model matrix
  y1 <- merged1$out1 ## select y-variable from NA-deleted dataset
  mod1 <- cv.glmnet(x1, y1, alpha = 1, family = "gaussian") ## run lasso model
  coef <- coef(mod1) ## get coef  
  
  list(merged1, x1, y1, mod1, coef)
  
}

# Clean Variables
## DVs
canada$dv_verylikely <- ifelse(is.na(canada$Q79)==F & canada$Q79==7,6,ifelse(is.na(canada$Q79)==F & canada$Q79<=6,7-canada$Q79,ifelse(is.na(canada$Q80)==F & canada$Q80==7,6,ifelse(is.na(canada$Q80)==F & canada$Q80<=6,7-canada$Q80,ifelse(is.na(canada$Q81)==F & canada$Q81==7,6,ifelse(is.na(canada$Q81)==F & canada$Q81<=6,7-canada$Q81,NA))))))
canada$dv_mostlikely <- ifelse(is.na(canada$Q79)==F & canada$Q79==7,7,ifelse(is.na(canada$Q79)==F & canada$Q79<=6,7-canada$Q79,ifelse(is.na(canada$Q80)==F & canada$Q80==7,7,ifelse(is.na(canada$Q80)==F & canada$Q80<=6,7-canada$Q80,ifelse(is.na(canada$Q81)==F & canada$Q81==7,7,ifelse(is.na(canada$Q81)==F & canada$Q81<=6,7-canada$Q81,NA))))))
canada$dv_missing <- ifelse(is.na(canada$Q79)==F & canada$Q79==7,NA,ifelse(is.na(canada$Q79)==F & canada$Q79<=6,7-canada$Q79,ifelse(is.na(canada$Q80)==F & canada$Q80==7,NA,ifelse(is.na(canada$Q80)==F & canada$Q80<=6,7-canada$Q80,ifelse(is.na(canada$Q81)==F & canada$Q81==7,NA,ifelse(is.na(canada$Q81)==F & canada$Q81<=6,7-canada$Q81,NA))))))
## Predictors
canada$region <- ifelse(canada$Q2>6,"West",ifelse(canada$Q2==6,"Ontario",ifelse(canada$Q2==5,"Quebec","Atlantic")))
canada$quebec <- ifelse(canada$region=="Quebec",1,0)
canada$west <- ifelse(canada$region=="West",1,0)
canada$agegroup <- ifelse(canada$Q3<35,"Age 18-34",ifelse(canada$Q3>=35 & canada$Q3<=44,"Age 35-44",ifelse(canada$Q3>=45 & canada$Q3<=54,"Age 45-54",ifelse(canada$Q3>=55 & canada$Q3<=64,"Age 55-64","Age 65+"))))
canada$age35 <- ifelse(canada$agegroup=="Age 35-44",1,0)
canada$age45 <- ifelse(canada$agegroup=="Age 45-54",1,0)
canada$age55 <- ifelse(canada$agegroup=="Age 55-64",1,0)
canada$age65 <- ifelse(canada$agegroup=="Age 65+",1,0)
canada$male <- ifelse(canada$Q4==1,1,0)
canada$university <- ifelse(canada$Q5>=9,1,0)
canada$married <- ifelse(canada$Q7==1,1,0)
canada$nonwhite <- ifelse(canada$Q11>1,1,0)
canada$frequentchurch <- ifelse(canada$Q13<=2,1,0)
canada$polinterest <- (5-canada$Q17)/4
canada$left <- ifelse(canada$Q18==1 | canada$Q18==3 | canada$Q18==4 | canada$Q18==5,1,0)
canada$left[canada$Q20==1 | canada$Q20==3 | canada$Q20==4 | canada$Q20==5] <- 1
canada$right <- ifelse(canada$Q18==2,1,0)
canada$right[canada$Q20==2] <- 1
canada$ideology1 <- canada$Q21
canada$approve_trudeau <- (4-canada$Q22)/3
canada$trudeau_therm <- canada$Q23_4/100
canada$vaccine_norms_guess <- canada$Q82/100
canada$underestimated_norms <- ifelse(canada$vaccine_norms_guess<0.61,1,0)
canada$overestimated_norms <- ifelse(canada$vaccine_norms_guess>0.79,1,0)
mean(canada$vaccine_norms_guess,na.rm=T)
quantile(canada$vaccine_norms_guess,c((1/3),(2/3)),na.rm=T)
canada$medium_tercile <- ifelse(canada$vaccine_norms_guess>=.68 & canada$vaccine_norms_guess<.81,1,0)
canada$high_tercile <- ifelse(canada$vaccine_norms_guess>=.81,1,0)
## Media Trust
cor.test(canada$Q15_1,canada$Q15_2) # r=.14
psych::alpha(cbind.data.frame(canada$Q15_1,canada$Q15_2)) ## alpha=.24
canada$total_media_trust <- ((4-canada$Q15_1)+(4-canada$Q15_2))/6
## Health Trust
psych::alpha(cbind.data.frame(canada$Q15_3,canada$Q15_5,canada$Q16_1,canada$Q16_2,canada$Q16_3,canada$Q16_8)) ## alpha=.91
canada$trust_healthgov <- ((4-canada$Q15_3)+(4-canada$Q15_5)+(4-canada$Q16_1)+(4-canada$Q16_2)+(4-canada$Q16_3)+(4-canada$Q16_8))/18
## Conspiracy Predispositions
psych::alpha(cbind.data.frame(canada$Q26_1,canada$Q26_2,canada$Q26_3,canada$Q26_4)) ## alpha=.85
canada$conspiracy <- ((5-canada$Q26_1)+(5-canada$Q26_2)+(5-canada$Q26_3)+(5-canada$Q26_4))/16
## Political Knowledge
canada$knowledge1 <- ifelse(canada$Q28==4,1,0)
canada$knowledge2 <- ifelse(canada$Q29==3,1,0)
canada$knowledge3 <- ifelse(canada$Q30==4,1,0)
canada$knowledge4 <- ifelse(canada$Q31==3,1,0)
canada$knowledge5 <- ifelse(canada$Q32==3,1,0)
psych::alpha(cbind.data.frame(canada$knowledge1,canada$knowledge2,canada$knowledge3,canada$knowledge4,canada$knowledge5)) ## alpha=.56
canada$knowledge <- (canada$knowledge1+canada$knowledge2+canada$knowledge3+canada$knowledge4+canada$knowledge5)/5
## Lagged DV
canada$lagdv_missing <- ifelse(canada$Q74==7,NA,7-canada$Q74)
canada$lagdv_verylikely <- ifelse(canada$Q74==7,6,7-canada$Q74)
canada$lagdv_mostlikely <- ifelse(canada$Q74==7,7,7-canada$Q74)
## Treatment
canada$norms_treatment <- factor(ifelse(is.na(canada$Q79)==F,"Control",ifelse(is.na(canada$Q80)==F,"Contention","Injunctive")),levels=c("Control","Contention","Injunctive"))
canada$injunctive_treatment <- ifelse(canada$norms_treatment=="Injunctive",1,0)
canada$contention_treatment <- ifelse(canada$norms_treatment=="Contention",1,0)

# Lasso Regression
df <- canada
lasso.dvmissing <- run_LASSO(df,dv_missing,lagdv_missing) # Health Trust, Pre
lasso.dvverylikely <- run_LASSO(df,dv_verylikely,lagdv_verylikely) # Health Trust, Pre
lasso.dvmostlikely <- run_LASSO(df,dv_mostlikely,lagdv_mostlikely) # Health Trust, Pre

# Table S3
table(canada$norms_treatment,canada$university)
chisq.test(canada$norms_treatment,canada$university)
canada$agegroup <- ifelse(canada$Q3<=34,"18 to 34",ifelse(canada$Q3>=35 & canada$Q3<=49,"35 to 49",ifelse(canada$Q3>=50 & canada$Q3<=64,"50 to 64","65 and older")))
table(canada$norms_treatment,canada$agegroup)
canada$age1834 <- ifelse(canada$agegroup=="18 to 34",1,0)
canada$age3549 <- ifelse(canada$agegroup=="35 to 49",1,0)
canada$age5064 <- ifelse(canada$agegroup=="50 to 64",1,0)
canada$age65 <- ifelse(canada$agegroup=="65 and older",1,0)
chisq.test(canada$norms_treatment,canada$age1834)
chisq.test(canada$norms_treatment,canada$age3549)
chisq.test(canada$norms_treatment,canada$age5064)
chisq.test(canada$norms_treatment,canada$age65)
table(canada$norms_treatment,canada$male)
chisq.test(canada$norms_treatment,canada$male)
table(canada$norms_treatment,canada$married)
chisq.test(canada$norms_treatment,canada$married)
table(canada$norms_treatment,canada$frequentchurch)
chisq.test(canada$norms_treatment,canada$frequentchurch)
table(canada$norms_treatment,canada$region)
canada$ontario <- ifelse(canada$region=="Ontario",1,0)
canada$atlantic <- ifelse(canada$region=="Atlantic",1,0)
chisq.test(canada$norms_treatment,canada$west)
chisq.test(canada$norms_treatment,canada$ontario)
chisq.test(canada$norms_treatment,canada$quebec)
chisq.test(canada$norms_treatment,canada$atlantic)
table(canada$norms_treatment,canada$right)
chisq.test(canada$norms_treatment,canada$hotspot)
table(canada$norms_treatment,canada$hotspot)
chisq.test(canada$norms_treatment,canada$right)
table(canada$norms_treatment,canada$left)
chisq.test(canada$norms_treatment,canada$left)
mean(subset(canada,norms_treatment=="Control")$ideology1,na.rm=T)
mean(subset(canada,norms_treatment=="Injunctive")$ideology1,na.rm=T)
mean(subset(canada,norms_treatment=="Contention")$ideology1,na.rm=T)
summary(aov(ideology1~norms_treatment,data=canada))
mean(subset(canada,norms_treatment=="Control")$knowledge*5,na.rm=T)
mean(subset(canada,norms_treatment=="Injunctive")$knowledge*5,na.rm=T)
mean(subset(canada,norms_treatment=="Contention")$knowledge*5,na.rm=T)
summary(aov(knowledge~norms_treatment,data=canada))
table(canada$norms_treatment,canada$nonwhite)
chisq.test(canada$norms_treatment,canada$nonwhite)
mean(subset(canada,norms_treatment=="Control")$trust_healthgov*3,na.rm=T)
mean(subset(canada,norms_treatment=="Injunctive")$trust_healthgov*3,na.rm=T)
mean(subset(canada,norms_treatment=="Contention")$trust_healthgov*3,na.rm=T)
summary(aov(trust_healthgov~norms_treatment,data=canada))
mean(subset(canada,norms_treatment=="Control")$total_media_trust*3,na.rm=T)
mean(subset(canada,norms_treatment=="Injunctive")$total_media_trust*3,na.rm=T)
mean(subset(canada,norms_treatment=="Contention")$total_media_trust*3,na.rm=T)
summary(aov(total_media_trust~norms_treatment,data=canada))
mean(subset(canada,norms_treatment=="Control")$lagdv_verylikely,na.rm=T)
mean(subset(canada,norms_treatment=="Injunctive")$lagdv_verylikely,na.rm=T)
mean(subset(canada,norms_treatment=="Contention")$lagdv_verylikely,na.rm=T)
summary(aov(lagdv_verylikely~norms_treatment,data=canada))
mean(subset(canada,norms_treatment=="Control")$lagdv_mostlikely,na.rm=T)
mean(subset(canada,norms_treatment=="Injunctive")$lagdv_mostlikely,na.rm=T)
mean(subset(canada,norms_treatment=="Contention")$lagdv_mostlikely,na.rm=T)
summary(aov(lagdv_mostlikely~norms_treatment,data=canada))
mean(subset(canada,norms_treatment=="Control")$lagdv_missing,na.rm=T)
mean(subset(canada,norms_treatment=="Injunctive")$lagdv_missing,na.rm=T)
mean(subset(canada,norms_treatment=="Contention")$lagdv_missing,na.rm=T)
summary(aov(lagdv_missing~norms_treatment,data=canada))

# Table S5/Figure 1
Reg.h1.verylikely.canada.robust <- lm_robust(dv_verylikely~injunctive_treatment+contention_treatment+trust_healthgov+lagdv_verylikely,data=canada)
ca1 <- effectsize(Reg.h1.verylikely.canada.robust)
Reg.h1.mostlikely.canada.robust <- lm_robust(dv_mostlikely~injunctive_treatment+contention_treatment+trust_healthgov+lagdv_mostlikely,data=canada)
Reg.h1.missing.canada.robust <- lm_robust(dv_missing~injunctive_treatment+contention_treatment+trust_healthgov+lagdv_missing,data=canada)

# Table S8/Figure 2
Reg.h1b.verylikely.canada.robust <- lm_robust(dv_verylikely~(injunctive_treatment*underestimated_norms)+(injunctive_treatment*overestimated_norms)+contention_treatment+trust_healthgov+lagdv_verylikely,data=canada)
Reg.h1b.canada.under <- summary(margins(Reg.h1b.verylikely.canada.robust,at=list(underestimated_norms=1)))
Reg.h1b.canada.over <- summary(margins(Reg.h1b.verylikely.canada.robust,at=list(overestimated_norms=1)))
Reg.h1b.mostlikely.canada.robust <- lm_robust(dv_mostlikely~(injunctive_treatment*underestimated_norms)+(injunctive_treatment*overestimated_norms)+contention_treatment+trust_healthgov+lagdv_mostlikely,data=canada)
Reg.h1b.missing.canada.robust <- lm_robust(dv_missing~(injunctive_treatment*underestimated_norms)+(injunctive_treatment*overestimated_norms)+contention_treatment+trust_healthgov+lagdv_missing,data=canada)

# Table S11
Reg.rq2.verylikely.canada.robust <- lm_robust(dv_verylikely~(injunctive_treatment*treat)+(contention_treatment*treat)+trust_healthgov+lagdv_verylikely,data=canada)
linearHypothesis(Reg.rq2.verylikely.canada.robust,c("treat=0","injunctive_treatment:treat=0","treat:contention_treatment=0"),test="F")
summary(margins(Reg.rq2.verylikely,at=list(treat=1)))
Reg.rq2.mostlikely.canada.robust <- lm_robust(dv_mostlikely~(injunctive_treatment*treat)+(contention_treatment*treat)+trust_healthgov+lagdv_mostlikely,data=canada)
linearHypothesis(Reg.rq2.mostlikely.canada.robust,c("treat=0","injunctive_treatment:treat=0","treat:contention_treatment=0"),test="F")
Reg.rq2.missing.canada.robust <- lm_robust(dv_missing~(injunctive_treatment*treat)+(contention_treatment*treat)+trust_healthgov+lagdv_missing,data=canada)
linearHypothesis(Reg.rq2.missing.canada.robust,c("treat=0","injunctive_treatment:treat=0","treat:contention_treatment=0"),test="F")

# Table S13
Reg.h1.verylikely.canada.robust.nocov <- lm_robust(dv_verylikely~injunctive_treatment+contention_treatment,data=canada)
Reg.h1.mostlikely.canada.robust.nocov <- lm_robust(dv_mostlikely~injunctive_treatment+contention_treatment,data=canada)
Reg.h1.missing.canada.robust.nocov <- lm_robust(dv_missing~injunctive_treatment+contention_treatment,data=canada)

# Table S16
Reg.h1b.verylikely.canada.robust <- lm_robust(dv_verylikely~(injunctive_treatment*underestimated_norms)+(injunctive_treatment*overestimated_norms),data=canada)
Reg.h1b.mostlikely.canada.robust <- lm_robust(dv_mostlikely~(injunctive_treatment*underestimated_norms)+(injunctive_treatment*overestimated_norms),data=canada)
Reg.h1b.missing.canada.robust <- lm_robust(dv_missing~(injunctive_treatment*underestimated_norms)+(injunctive_treatment*overestimated_norms),data=canada)

# Table S19
Reg.rq2.nocov.verylikely.canada.robust <- lm_robust(dv_verylikely~(injunctive_treatment*treat)+(contention_treatment*treat),data=canada)
linearHypothesis(Reg.rq2.nocov.verylikely.canada.robust,c("treat=0","injunctive_treatment:treat=0","treat:contention_treatment=0"),test="F")
Reg.rq2.nocov.mostlikely.canada.robust <- lm_robust(dv_mostlikely~(injunctive_treatment*treat)+(contention_treatment*treat),data=canada)
linearHypothesis(Reg.rq2.nocov.mostlikely.canada.robust,c("treat=0","injunctive_treatment:treat=0","treat:contention_treatment=0"),test="F")
Reg.rq2.nocov.missing.canada.robust <- lm_robust(dv_missing~(injunctive_treatment*treat)+(contention_treatment*treat),data=canada)
linearHypothesis(Reg.rq2.nocov.missing.canada.robust,c("treat=0","injunctive_treatment:treat=0","treat:contention_treatment=0"),test="F")

# Figure S1 (Canada)
Reg.h1b.verylikely.canada.flex <- interflex("kernel",data=canada,Y="dv_verylikely",D="injunctive_treatment",X="vaccine_norms_guess",Z=c("contention_treatment","trust_healthgov","lagdv_verylikely"),na.rm=T,CI=T)